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Abstract 

For the Nonlinear Shrodinger Equation with disorder it was found nu- 
merically that in some regime of the parameters Anderson localization is 
destroyed and subdiffusion takes place for a long time interval. It was 
argued that the nonlinear term acts as random noise. In the present 
work the properties of this effective noise are studied numerically. Some 
assumptions made in earlier work were verified, and fine details were ob- 
tained. The dependence of various quantities on the localization length of 
the linear problem were computed. A scenario for the possible breakdown 
of the theory for a very long time is outlined. 

1 Introduction 

The Nonlinear Schrodinger Equation (NLSE) [1J in a random potential takes 
the form of 

id t ^ = H^ + p\^\ 2 ^, (i) 

where Hq is the linear part with a disordered potential, which on a lattice takes 
the form of 

H ip(x) = - (tl>(x + 1) + if>(x - 1)) + e(x)i/>(x). (2) 

In this work it is assumed that e (x) are identical independent random variables 
(i.i.d) uniformly distributed in the interval of [— . 

The NLSE was derived for a variety of physical systems under some ap- 
proximations. It was derived in classical optics where tp is the electric field by 
expanding the index of refraction in powers of the electric field keeping only the 
leading nonlinear term [2]. For Bose-Einstein Condensates (BEC), the NLSE is 
a mean field approximation where the term proportional to the density ap- 
proximates the interaction between the atoms. In this field the NLSE is known 
as the Gross-Pitaevskii Equation (GPE) O HI El [7] . It is well known that 
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in ID in the presence of a random potential with probability one all the states 
are exponentially localized [H [9l [10] . Consequently, diffusion is suppressed and 
in particular a wavepacket that is initially localized will not spread to infinity. 
This is the phenomenon of Anderson localization [11J . The problem defined 
by (TT]) is relevant for experiments in nonlinear optics, for example disordered 
photonic lattices |12[ [T5] , where Anderson localization was found in presence of 
nonlinear effects as well as experiments on BECs in disordered optical lattices 
[II [T^ [TJ [T71 [II [II [201 13 111 123] • The interplay between disorder and nonlin- 
ear effects leads to new interesting physics [101 [HI [Ml 12S1 1211 ■ In spite of the 
extensive research, many fundamental problems are still open (see recent review 
to be published in Nonlinearity [28J). In particular there is disagreement be- 
tween the analytical and the numerical results [2lfl |501 15 3 1521 1551 154 1 15 ^ 156 ] 157] . 

A natural question is whether a wave packet that is initially localized in 
space will indefinitely spread for dynamics controlled by (JTJ. A simple argu- 
ment indicates that spreading will be suppressed by randomness. If unlimited 
spreading takes place the amplitude of the wave function will decay since the 
l 2 norm is conserved. Consequently, the nonlinear term will become negligi- 
ble and Anderson localization will take place as a result of the randomness as 
conjectured by Frohlich et al [55]. Contrary to this intuition, based on the 
smallness of the nonlinear term resulting from the spread of the wave function, 
it is claimed that for the kicked-rotor a nonlinear term leads to delocalization 
if it is strong enough [5J5] • It is also argued that the same mechanism results in 
delocalization for the model ([I]) with sufficiently large j3, while, for weak nonlin- 
earity, localization takes place [39, 40] . Recently, it was rigorously shown that 
the initial wavepacket cannot spread so that its amplitude vanishes at infinite 
time, for large enough /? [41]. It does not contradict spreading of a fraction 
of the wavefunction. Indeed, subdiffusion was found in numerical experiments 
[391 EtOl Et2l Et5I EH] . It was also argued that nonlinearity may enhance discrete 
breathers [26, 27J. In conclusion, it is not clear what is the long time behavior 
of a wave packet that is initially localized, if both nonlinearity and disorder are 
present |28| . The major difficulty in numerical resolution of this question is in- 
tegration of (TT]) to large time. Most researchers who run numerical simulations 
use a split-step method for integration, however it is impossible to achieve con- 
vergence for large times, and therefore some heuristic arguments assuming that 
the numerical errors do not affect the results qualitatively, are utilized [39"1 143], 
Moreover the problem is chaotic, therefore the trajectories that are found are 
not the actual trajectories and it is argued that it does not affect the statistical 
results. 

Recent rigorous arguments |291 130| in the limit of strong disorder combined 
with perturbation theory [53 [551 EES] indicate that it is unlikely that sub-diffusion 
persists forever and the asymptotic growth is at most logarithmic in time. Also 
other recent work based on a scaling theory [55] and phase space considerations 
[361 0S] lead to similar indications. It is clear that there is a substantial regime 
in time and parameters where sub-diffusion may hold and the purpose of the 
present work is to analyze the dynamics in this regime. 

Our analysis based on |43] 144] . is conveniently expressed expanding the wave- 
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function 

i/j(x,t) = Y J c n {t)u n {x)e- iE - t (3) 

n 

where u n are the eigenfunctions of H$ typically falling off exponentially: 

e -\x n -x\/£ 

u„(x) « — <p{x) (4) 

where p(x) is a random function of order unity. The localization center is x n . 
The c n (t) satisfy 



id t c n {t) = /3 V^^e^+^-^-^c^c^c^ = F n (t) (5) 

and 



mi ,7ri2 ,m3 



C 1,m2,m3 =J2 U n( X ) U m 1 (x)u m2 ( X )u m3 {x). (6) 

x 

In [43l HI] it is argued that F n (t) behaves as random noise with rapidly decay- 
ing correlation functions. The implications are analyzed in Sec. 2 and tested 
numerically in Sec. 3. A scenario for the breakdown of the effective noise the- 
ory is outlined in Sec. 4. The results are summarized and open question are 
presented in Sec. 5. 



2 The effective noise theory 

The theory of SKFF (Skokos, Krimer, Komineas and Flach [331 EH) assumes 
for spreading to take place to the region where the n-th state is localized from 
the region where the states mi, 11J2, "13 have a large amplitude: 

|c mi | 2 w |c m2 | 2 w |c ma | 2 «p (7) 

while 

\Cn\ 2 < P (8) 

It is assumed that the RHS of ([5]) is a random function denoted by F n (t). We 
turn to estimate its typical behavior. First we note that the overlap sums ([6]) 
are random functions. Within the scaling theory for localization one expects 
that for sufficiently weak disorder their various moments are determined by the 
localization length. For the case where all indices (n, mi, m2, TO3) are identical 
the average is just the inverse participation ratio what is proportional to For 
the general case the scaling theory suggests it is a function only of £. Experience 
with scaling theories leads us to assume it is a power of £. Therefore we try the 
form, 

^ V rmi,m a> m 8 ) = C r(l)^-»» j (9) 
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and for the second moment we try to fit to, 

< |y™i^2,m3|2 >= C r(2)^-2i» < ( 1Q ) 

Here Cq and are constants and < .. > is an average over realizations. 
We note that when the m, and n are all different the average of the overlap 
integrals vanishes. We should note that the localization length £ is actually 
energy dependent. For weak disorder in the center of the band, £ ~ W~ 2 [4"Tll48| . 
this relation holds for most energies in the energy band [47 . In what follows we 
will estimate the values of r\\ and 772 for various disorder strengths and for various 
sites (x n ,x mi ,x m2 ,x m3 ), which are within the localization length. Otherwise 
the sum ^ is negligible. It is not obvious that both © and (fTU)) will scale 
in this way although it is expected from the scaling theory of localization, that 
this is the case for sufficiently weak disorder, namely large £. We demonstrare 
that this is indeed the case and there is a typical magnitude of the value of the 
of the overlap sum (O and it scales as, 

V = CiC* (11) 

where C\ is a constant. Here and in what follows we denote by £ the localization 
length in the center of the band. 

Making the assumption that F n is random combined with ([7]) , the sum 
on the RHS of ([5]) consists of the order of £ 3 terms, at least for weak disorder. 
These are rapidly oscillating in time, and it is a nonlinear function of the c mi (t). 
Therefore it is suggestive that it can be considered random. This assumption 
will be tested in detail in subsection 13. II The RHS of © is assumed to take the 
form [H] 

F n = VVPp^Ut) = ^VPp^Ut) (12) 
where C± is a constant and 

V = A j3^ a P (13) 
is proportional to the number of "resonant modes", namely ones that strongly 
affect the dynamics of the state n. Although it is reasonable to assume that the 
number of resonant modes is proportional to the density p a strong argument for 
it is missing, nevertheless it is consistent with all numerical results |44l 143) . We 
assume here the form (fT5]) where A$ is a constant independent of j3 and £. In the 
end of this section we argue that within these assumption 7 = 1 in agreement 
with the assumption of [HI S3] ■ The value of a is estimated numerically (see 
subsection 13. 3|) . Under these assumptions ([5]) reduces to: 

id t c n {t) = F n (t) (14) 

Assuming F n (t) can be considered random with rapidly decaying correlations, 
in particular we assume that the distribution function of /„ (t) is stationary 
and the integral of correlation function C (t') = (f (0) f (t')), where (..) is the 
average over the random potential, converges. Integration results in 

c n {t)=-i'^VfSp s ' 2 J*dt'f n (f) (15) 
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Integrating over a time interval which is sufficiently large yields: 



< M*)| 2 >= |^Wi = A.Al^^p^-^t (16) 

where A\ is a constant. The value of < |c„(i)| 2 > increases with time and 
equilibrium is achieved when it takes the value p. Transitions between states 
of the type of n (states with small amplitude) are ignored in this model. The 
required time for equilibration is 

r =*FV (17) 

where we define 

B = A.Al^+^e^ 2 ^ 2 (18) 

The equilibration time T varies slowly compared to t (see discussion after ([24)) ). 
In other words there is a separation of time scales. On the time scale T the 
system seems to reach equilibrium by a diffusion process and the density becomes 
constant in a region that includes the site n. Hence on this time scale it seems 
to equilibrate. On a longer time scales, there is an even longer equilibration 
time scale, and the resulting diffusion is even weaker. The consistency of the 
argument results of the fact that ^ — > for t — > oo. Therefore it is assumed 
that the variations of p and T are slow on the scale of t . This assumption is 
checked in the end of this section. The resulting diffusion coefficient is 

D = C^ = CBp 4 (19) 

where C is a constant. The assumption is that the nonlinear term generates a 
random walk with the characteristic steps T and £ in time and space. At time 
scales t~^> T, there is diffusion and 

M 2 = Dt, (20) 

where Mi = ^ x \ip (x, t) | 2 and the variance M 2 = ^ (x — Mi) 2 \ip (x, t) | 2 are 
the first and second moments. Since the second moment M 2 is inversely pro- 
portional to p 2 one finds 

\ = A 2 CB P H (21) 
P 

where A 2 is a constant. Therefore 



P 

The second moment satisfies: 



4 = {A 2 CBt) 1/3 . (22) 



M 2 = -L-(CBt)V 3 (23) 



AT 
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and 



^ _ i _ c 2 / 3 A 2 2 /3 et 2 / 3 _ ce 

The density p and the equilibration time T change with time as p ~ and 
T ~ t§. Therefore for ^ ^ t'i and ^ ~ <T§. First note that in the long 
time limit t — > oo both derivatives vanish and -C ^jjr. Therefore for the 
derivation of the equilibration time /? can considered constant and on long scales 
of spreading T and D can be considered constant. Therefore the theory is 
consistent for large t. Since in the NLSE j3 appears only via the combination 
(3\ip(x)\ 2 , it can appear in (fT8|) and (fT9]) only in the power 4 (that is in the 
combination /3 4 p 4 ) therefore 7=1. 

In the next section this theory will be tested numerically. 



3 Numerical tests for the effective noise theory 

In this section the theory presented in section 2 is tested numerically. In sub- 
section 13.11 the distribution of the F n (t) is computed, in subsection 13.21 the 
first moments of the overlap sums are calculated while in subsection 13.31 the 
dependence of the second moment M 2 of (f2"3"|) on £ is evaluated. 



3.1 Statistical properties of F n (t) 

In this subsection the statistical distribution of F n (t) is explored. For this 
purpose the time dependent NLSE ([1} was solved numerically for a finite lattice 
of N sites, for Nr realizations of the random potential e (x) and for W — 4. The 
wavefunction ip (x, t) at time t was calculated for a single site excitation namely 
the initial condition ip (x, 0) = S Xy o using the split step method. The details of 
the numerical calculation are presented in the appendix. The expansion ([3]) of 
ip in terms of eigenfunctions of the linear problem @ yields, 

id t c n (t)=J2P\ip(x,t)\ 2 ip(x,t)u n (x)e ltE ™ = F n {t). (25) 

X 

This equation was used to calculate F n (t) numerically for a lattice of N sites. 
In order to check whether F n (t) can be considered as noise we calculated its 
power spectrum and auto-correlation function. First we present results obtained 
for times up to t = 10 5 for /3 = 1, W = 4 (f « 6.4), N = 1024 for a single site 
excitation at t = 0. The calculation was preformed for Nr — 50 realizations. 
For nearly all these realizations it was found that the second moment grows as 
M 2 oc i 1 / 3 in agreement with the results of |40[ 231 S3] ■ We focus first on such 
realizations and present the results for a specific realization in Fig. [1] 
The power spectrum is 

S n H = F n (w) 2 , (26) 
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where 



F n (w) = Jim -= \ F n (t) ■ e(~ iu)t Ut. 




(27) 



o 



It is plotted for some realization in Fig. [TJa for n=0. It exhibits a peak around 
| wo | ~ 1-72 and its width is Aco m 0.1. The finite width is characteristic of 
noise. Also the Fourier transform of 



will exhibit a wide power spectrum near u = 0, with the width of Au; that is 
characteristic of noise. The auto-correlation function of F n (t) is 



For F n (t) we define the auto-correlation function C n (r) that is just (|29l) 
with F n (t) replaced by F n (t) . In Fig. Q] .b we plot d R) = Re (C„ (r)) for 
n = while in Figfl] c the zoomed version is plotted. Note an oscillation of 
frequency of the order \ljq\ ~ 1.72 that is superimposed on the function. In the 

corresponding plots of C n R ^ = Re [c n (r)j , presented in Fig[TJd and Fig[TJe, 
one does not find this oscillation. Behavior of the imaginary part of the auto- 
correlation function CrP — Im (c n (rfj is similar (see FigfTJf). All results 
presented in FigfT] are for n = 0. Similar results were found also for n = 3 
and n = 15. We see that the auto-correlation function decays by 2 orders of 
magnitude on the scale of At ss 140 (of the order of 2ir/ Auj ~ 65). Therefore the 
correlation of F n (t) behaves as the one of noise with short time correlations. 
For realizations where the growth of the second moment M% ~ i 1 / 3 was not 
found, the power spectrum was found to be substantially narrower by 2 orders 
of magnitude. The calculations were repeated for f3 = 2 where similar results 
were found, and for j3 = 0.5. For the latter case the number of realizations where 
it was found that the second moment grows like t 1 ' 3 is substantially smaller 
than for /3 = 1 or j3 = 2. In all cases where the width of the power spectrum 
was small the typical growth of the secomd moment M% ~ i 1 / 3 was not found 
and vice versa. This demonstrates the strong relation between the effective noise 
behavior and the diffusive growth of the second moment. It also demonstrates 
the different behavior of various realizations of the randomness. 

We turn now to test the distribution of F n (t). For this purpose we sample 
F n (t) for a sequence of points separated by t a > At , that is for points where 
the values of F n (t) are uncorrelated, and compute the distribution of F n (k ■ t a ) 
for k = (1, 2, ..K). The results are presented in Fig. Hfor t = 10 5 , t a = 200, 
K = 500. 



F n (t) = F n (t) • e 



(28) 



C n (r) = F n (t)-F*(t + T) 



(29) 



where bar denotes time average g (t) = lim^^ i J* g (t) dt. 
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Figure 1: The corelation C n (t) and power spectrum S n (to) of F n (t) for W = 4 
, p = 1, N = 1024, t = 10 5 ,n = 0. (a) The Power Spectrum S (u>), (b) 

( R) ( R) 

The auto-correlation function Cq (t). (c) The zoomed Cq (r), (d) The 

auto-correlation function Cq (r), (e) the zoomed Cq (t). (f) the zoomed 

cjp (r)[see text]. 





Figure 2: The distribution of Y = F,[ R) (k ■ t a ) where k = (1, 2, ..K) , K = 500, 
t a = 200 , t — 10 5 and the bin size 0.0596 . (a) For the same realization used in 
Fig. 1 . (b) The distribution of values found for all Nr — 50 realizations. 



3.2 Estimate of scaling of the matrix elements v^" 1 " 12 " 13 
with £ 

The overlap sum V™ 1 '™ 2 '™" 3 is a random function. In this subsection the scaling 
of its typical values with the maximal localization length [47J 
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(30) 



is evaluated. This relation holds in the limit of weak disorder. In the numerical 
calculations presented in this paper we vary W as the control parameter and the 
localization length is calculated from (|30|) . The estimate (|30|) is a reasonable ap- 
proximation for W < 5.5 or £ > 3.15 as was checked explicitly (and used) in this 
subsection. We note that the v™ 1,m2,ni3 take values of substantial magnitude 
when all the centers of localization of the states u n , u mi , u m2 , u m3 are within a 
distance £. Only such overlap sums are considered. The average of the overlap 
sums over realizations vanishes unless (n,mi, m2,rn^) consists of two pairs of 
identical values , n = mj and 7712 = m 3 and all permutations. We calculated 

/|^i,m 2 ,m 3 |2^ and (ymx,m*,m 3 

realizations) while x 
are varied. Assuming (V™ 1 '™ 2 '™ 3 ) ~ and /|V^i."»a,m 3 | 2 \ _ £-?V2 while 

the variance (^(y™ l < m2 < m3 ) 2 \ - (V^^^,™^ 2 scales as ^-2» )3; we estimate these 

exponents from Figures like Fig. [3] . We conclude that 771 « r\ 2 ~ ?73 ~ 1 • 
Therefore the typical magnitude of the random variable y,™ 1 '™ 2 '™ 3 scales as 
(fTTj) with rj = 1. Although this result is expected from the scaling theory of 
localization, it is not obvious appriory. In particular it is not clear what is the 
effect of cancellations of various terms resulting of opposite signs. 

For £ <C 11 we could not obtain smooth curves of V™ 1 '™ 2 '™ 3 . The reason is 



(where (•) denotes average over Nr — 5000 
mi , «u m2 , x m3 are fixed fractions of £, while £ (and W) 
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(g) 



x=ln© 
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Figure 3: A log- log plot of (b) y = In I V ' 3 ' 3 V (r) y 

2\ / , , \ 2\ 
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and (g) y = In I ( yV ' 3 ' 3 J y — y Vjj ' 3 ' 3 y J as a function of .t = /ri(£) , for 

the parameters iV = 512 , Nr = 5000. The localization length varies in the 
interval 11 < £ < 103. The least square fit leads to r\\ — 1.039 , r\2 — 0.958 and 
773 = 0.853 respectively . 



that the centers of localization equal to the integer part of i/a where a 

is fixed and £ varies. For small £ the jumps in V^ ll ' rn:>,m3 are significant, since 
£ does not cover many integers. The results obtained indicate that scaling of 
the overlap sums as £ _1 holds also for values £ < 11. In summary for a crude 
evaluation one can assume (|11[) holds with 77 = 1. 



3.3 The scaling of the second moment M 2 with £ ( and /?) 

In this subsection we will estimate the exponent a defined in (j!3[) . For this 
purpose we write (|23p in the form 

M 2 = At* (31) 

with 

A = A A C (32) 

where v — | (a — r\ + 1) (see (fT8)) ) while is a constant independent of £. We 
used the split step method to obtain ip (x,t) for different realizations (Nr = 30) 
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-1 1 1 ' ' ' ' 1 

0.5 1 1.5 2 2.5 3 3.5 

x-lng) 

Figure 4: The dependence of ^4 defined by (|31l) and (|32p for /3 = 1 (blue circles) 
and for (3 — 3 (red squares) on £. We denote y = In (A) and a; = Zn(£) . From 
the least square fit we find v — 1.684 for (3 — 1 (blue) and v — 1.395 for j3 = 3 
(red). 

and computed V> until i = 10 . Only realizations which satisfied M2 ~ at 
some stage of the calculation were taken into account. This was the case for 
nearly all the Nr realizations for £ > 7 and j3 < 4. In the other regimes it was 
not satisfied for a significant number of realizations. Fixing /3 we estimate v 
from plots like Fig. HI For 1 < f3 < 3.5 using the fact that n « 1 we find that 
for 1.235 < v < 1.71 for various values of /3. The exponent a of ([TBI takes the 
values 1.85 < a < 2.56 . We note the strong uncertainty of v and a. These 
results indicate that It is an estimate of the order of magnitude but not 

a verification of this power law. 

4 Possibility for the breakdown of the effective 
noise theory 

For the effective noise theory it is essential that F n (t) can be considered random. 
For this the number of terms in the sum ([5]) that resonate with n should be 
large, namely V should not be too small. The density p and therefore V decrease 
with time. If V is very small there may be a situation that as a result of 
fluctuations, the sum ([5]) is dominated just by one term and therefore it is 
effectively quasi periodic. If spreading is a result of the randomness of F n , it 
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will stop then. Let us first estimate the time scale required to spread so that 
V s» 1. For this purpose let us write (fT3|) in the form 



V « Ae a /J (33) 

where A = Aq/3. Since p decreases with time t there is a time scale when V will 
become very small. Assuming the constants are of the order of unity, using (|18|) 
and ([21]) the time t* when V ~\ satisfies 

e a r « 1 (34) 

r£2(<*-77+l)£*l 3 

or 

f (l«*+f „ (35) 

resulting in 

f « ^( 4 «+ 2 ("- 1 )) (36) 

for 1.85 < a < 2.56 and rj = 1 

t* « £ 5 (37) 

where 7.4 < <5 < 10.24 

The time required for, V -C 1, when the effective noise theory may fail is even 
larger. 



5 Summary and conclusions 

The effective noise theory was introduced in [39J and was further developed 
in |431 14171 [U]. It was found to be consistent with the numerical results in 
some regimes. In Section [5] our interpretation of this theory was presented. 
In section [3] the details of this theory were tested numerically. In particular 
the distribution of the effective driving F n defined in ([5]) was studied . The 
correlation function was calculated as well and was found to be characterized 
by a wide power spectrum and rapid decay with time. These were found only for 
realizations where subdiffusion with the second moment growing as t 1 / 3, is found, 
indicating the relation between this spreading and the approximation of F n as 
effective noise. These results are purely numerical and support the effective 
noise theory. An obvious challenge is to obtain these results analytically. We 
determined that the behavior A w (see (f3"2")) ). with 1.235 < v < 1.71 is a 
reasonable approximation. From this we conclude that the dependence of V on 
£ (fT5|) is controlled by the exponent 1.85 < a < 2.56. Although £ varied over 
one decade and the evaluation of the exponent is crude we believe it may give 
the correct order of magnitude. 

We turn to speculate how the effective noise theory may break down for a 
long time scale. Assuming the effective noise theory holds for long time, V of 
(|13p becomes extremely small, consequently the number of terms in the sum 
([5]) that contribute significantly may become of order unity and F n may turn to 
be quasi periodic rather that random. Therefore there is a time scale t* given 
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by the estimate (f3"51 so that for t > t* the effective noise theory is invalid. For 
such long time a sequence of peaks may replace the continuous region of the 
power spectrum in Fig l.a. If localization is destroyed by the effective noise 
F n , it is reasonable to expect localization or spreading slower than subdiffusion 
(say logarithnic in time) on time scale t* and larger. Existence of such a time 
scale is consistent with |29[ [BUI I3"T1 |3"5I [551 143] , The scaling arguments used 
here should improve when the localization length £ becomes large but then t* 
becomes extremely large and it is impossible to explore numerically the scenario 
for the breakdown of the effective noise theory outlined in Sec. 0] Such a 
scenario may enable to reconcile the numerical results where subdiffusion is 
found [2SJ S3 Ell EH H3J H5 with the analytical results predicting asymptotic 
spreading that is at most logarithmic [25J [251 1301 131] • These points should be 
subject of future research. 
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Appendix: some details of the numerical calcula- 
tions 

We used the split step method to obtain the time evolution starting from the 
initial wavefunction. The lattice size N used is 512 or 1024. The reason we 
used the relativity large lattice is because we wanted to avoid boundary effects, 
namely we required the wavefnction amplitude to be smaller than 10~ 12 on 
the boundary. The time step used in the split step method is dt = 0.1. We 
used this time step because it is small enough relative to the time scales in the 
system at hand and large enough in order to complete the numerical calculation 
in reasonable time. It is the smallest time step used in |431 144j. The initial 
condition used is a single site excitation in the middle of the lattice denoted by 
x n = namely, tji (x,t = 0) = 5 X fl. 
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